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Abstract 

Heterogeneous materials abound in nature and man-made situations. Examples include porous 
media, biological materials, and composite materials. Diverse and interesting properties exhibited 
by these materials result from their complex microstructures, which also make it difficult to model 
the materials. Yeong and Torquato [Phys. Rev. E 57, 495 (1998)] introduced a stochastic op- 
timization technique that enables one to generate realizations of heterogeneous materials from a 
prescribed set of correlation functions. In this first part of a series of two papers, we collect the 
known necessary conditions on the standard two-point correlation function 5*2 (r) and formulate a 
new conjecture. In particular, we argue that given a complete two-point correlation function space, 
52 (r) of any statistically homogeneous material can be expressed through a map on a selected set of 
bases of the function space. We provide new examples of realizable two-point correlation functions 
and suggest a set of analytical basis functions. We also discuss an exact mathematical formulation 
of the (re) construction problem and prove that 52(r) cannot completely specify a two-phase het- 
erogeneous material alone. Moreover, we devise an efficient and isotropy-preserving construction 
algorithm, namely, the Lattice- Point algorithm to generate realizations of materials from their two- 
point correlation functions based on the Yeong- Torquato technique. Subsequent analysis can be 
performed on the generated images to obtain desired macroscopic properties. These developments 
are integrated here into a general scheme that enables one to model and categorize heterogeneous 
materials via two-point correlation functions. We will mainly focus on the basic principles in this 
paper. The algorithmic details and applications of the general scheme are given in the second part 
of this series of two papers. 

PACS numbers: 05.20.-y, 61.43.-j 
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I. INTRODUCTION 



A heterogeneous material (medium) is one that composed of domains of different materi- 
als or phases (e.g., a composite) or the same material in different states (e.g., a polycrystal). 
Such materials are ubiquitous; examples include sandstones, granular media, animal and 
plant tissue, gels, foams and concrete. The microstructures of heterogeneous materials can 
be only characterized statistically via various types of n-point correlation functions [1]. The 
effective transport, mechanical, and electromagnetic properties of heterogeneous materials 
are known to be dependent on an infinite set of correlation functions that statistically char- 
acterize the microstructure 

Reconstruction of heterogeneous materials from a knowledge of limited microstructural 
information (a set of lower-order correlation functions) is an intriguing inverse problem 

mm 

. An effective reconstruction procedure enables one to generate accurate structures 
and subsequent analysis can be performed on the image to obtain macroscopic properties 
of the materials; see, e.g., Ref. Gj. This provides a nondestructive means of estimating 
the macroscopic properties: a problem of important technological relevance. Another useful 
application is reconstruction of a three-dimensional structure of the heterogeneous material 
using information extracted from two-dimensional plane cuts through the material Such 
reconstructions are of great value in a wide variety of fields, including petroleum engineer- 
ing, biology and medicine, because in many cases one only has two-dimensional information 
such as a micrograph or image. Generating realizations of heterogeneous materials from 
a set of hypothetical correlation functions is often referred to as a construction problem. 
A successful means of construction enables one to identify and categorize materials based 
on their correlation functions. One can also determine how much information is contained 
in the correlation functions and test realizability of various types of hypothetical correla- 
tion functions. Furthermore, an effective (re) construction procedure can be employed to 
investigate any physical phenomena where the understanding of spatiotemporal patterns is 



fundamental, such as in turbulence 



A popular (re) construction procedure is based on the use of Gaussian random fields: 
successively passing a normalized uncorrelated random Gaussian field through a linear and 
then a nonlinear filter to yield the discrete values representing the phases of the structure. 
The mathematical background used in the statistical topography of Gaussian random fields 



3 



technique 



Torquato a: 



was originally established in the work of Rice [8| , l9 | . Many variations of this method have 
been developed and applied since then 10, III, 12, 13]. The Gaussian-field approach assumes 
that the spatial statistics of a two-phase random medium can be completely described by 
specifying only the volume fraction and standard two-point correlation function 5*2 (r), which 

B'ves the probability of finding two points separated by vector distance r in one of the phases 
]. However, to reproduce Gaussian statistics it is not enough to impose conditions on the 
first two cumulants only, but also to simultaneously ensure that higher-order cumulants 
vanish |l^ . In addition, the method is not suitable for extension to non-Gaussian statistics, 
and hence is model dependent. 

Recently, Torquato and coworkers have introduced another stochastic (re) construction 

In this method, one starts with a given, arbitrarily chosen, initial 
configuration of random medium and a set of target functions. The medium can be a dis- 
pe.io„ of pa..,e-li.e build., block. ,15] o. .o.e .en.^1. a „ ..age g fl Q. 
The target functions describe the desirable statistical properties of the medium of inter- 
est, which can be various correlation functions taken either from experiments or theoretical 
considerations. The method proceeds to find a realization (configuration) in which calcu- 
lated correlation functions best match the target functions. This is achieved by minimizing 
the sum of squared differences between the calculated and target functions via stochastic 
optimization techniques, such as the simulated annealing method [l^. This method is appli- 
cable to multidimensional and multiphase media, and is highly flexible to include any type 
and number of correlation functions as microstructural information. It is both a general- 
ization and simpliflcation of the aforementioned Gaussian-fleld (re) construction technique. 
Moreover, it does not depend on any particular statistics [1,|. 

There are many different types of statistical descriptors that can be chosen as target 



functions [Ij; the most basic one is the aforementioned^ two-point correlation function 52 (r), 
which is obtainable from small-angle X-ray scattering [17]. However, not everyhypothetical 
two-point correlation function corresponds to a realizable two-phase medium [l| . Therefore, 
it is of great fundamental and practical importance to determine the necessary conditions 



that realizable two-point correlation functions must possess |18l. |19|. Shepp showed that con- 
vex combinations and products of two scaled autocovariance functions of one-dimensional 
media (equivalent to two-point correlation functions; see deflnition below) satisfy all known 
necessary conditions for a realizable scaled autocovariance function 20| . More generally, we 
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will see that a hypothetical function obtained by a particular combination of a set of realiz- 
able scaled autocovariance functions corresponding to dimensional media is also realizable. 

In this paper, we generalize Shepp's work and argue that given a complete two-point 
correlation function space, 5*2 (r) of any statistically homogeneous material can be expressed 
through a map on a selected set of bases of the function space. We collect all known necessary 
conditions of realizable two-point correlation functions and formulate a new conjecture. We 
also provide new examples of realizable two-point correlation functions and suggest a set 
of analytical basis functions. We further discuss an exact mathematical formulation of 
the (re) construction problem and show that 5'2(r) cannot completely specify a two-phase 
heterogeneous material alone, apart from the issue of chirality. Moreover, we devise an 
efficient and isotropy-preserving construction algorithm to generate realizations of materials 
from their two-point correlation functions. Subsequent analysis can be performed on the 
generated images to estimate desired macroscopic properties that depend on 5'2(r), including 



both linear 
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261] and nonlinear 27|, 1281] behavior. These developments are 



integrated here into a general scheme that enables one to model and categorize heterogeneous 
materials via two-point correlation functions. Although the general scheme is applicable in 
any space dimension, we will mainly focus on two-dimensional media here. In the second 
part of this series of two papers 29|, we will provide algorithmic details and applications of 
our general scheme. 

The rest of this paper is organized as follows: In Sec. II, we briefly introduce the basic 
quantities used in the description of two-phase random media. In Sec. Ill, we gather all 
the known necessary conditions for realizable two-point correlation functions and make a 
conjecture on a new possible necessary condition based on simulation results. In Sec. IV, we 
propose a general form through which the scaled autocovariance functions can be expressed 
by a set of chosen basis functions and discuss the choice of basis functions. In Sec. V, 
we formulate the (re) construction problem using rigorous mathematics and show that 5*2 (r) 
alone cannot completely specify a two-phase random medium. Thus, it is natural to solve the 
problem by stochastic optimization method (i.e., simulated annealing). The optimization 
procedure and the Lattice-Point algorithm are also discussed. In Sec. VI, we provide several 
illustrative examples. In Sec. VII, we make concluding remarks. 
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II. DEFINITIONS OF n-POINT CORRELATION FUNCTIONS 



The ensuing discussion leading to the definitions of the n-point correlation functions 
follows closely the one given by Torquato Consider a realization of a two-phase random 
heterogeneous material within d-dimensional Euclidean space R'^. To characterize this binary 
system, in which each phase has volume fraction 0j [i = 1, 2), it is customary to introduce 
the indicator function /('^(x) defined as 

(1) 

[ 0, X e Vi, 

where G R'' is the region occupied by phase i and Vi G M*^ is the region occupied by the 
other phase. The statistical characterization of the spatial variations of the binary system 
involves the calculation of n-point correlation functions: 

^«(Xi,X2, • • • ,XJ = (/«(X1)/«(X2) ■ --I^'K^n)) , (2) 

where the angular brackets (■ ■ ■) denote ensemble averaging over independent realizations of 
the random medium. 

For statistically homogeneous media, the ra-point correlation function depends not on the 
absolute positions but on their relative displacements, i.e., 

5«(xi,X2, ■ ■ ■ ,xO = 5»(xi2, ■ ■ ■ ,xiO, (3) 

for all n > 1, where Xjj = Xj — Xj. Thus, there is no preferred origin in the system, which in 
Eq. ([3]) we have chosen to be the point Xi. In particular, the one-point correlation function 
is a constant everywhere, namely, the volume fraction (pi of phase i, i.e., 

^P = (/«(x)> = 0„ (4) 

and it is the probability that a randomly chosen point in the medium belongs to phase i. For 
statistically isotropic media, the n-point correlation function is invariant under rigid-body 

(i) 

rotation of the spatial coordinates. For n < d, this implies that Sn depends only on the 
distances Xij = |xjj| (1 < i < j < n). For > ci-l- 1, it is generally necessary to retain vector 
variables because of chirality of the medium. 
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The two-point correlation function S'2*^(xi,X2) defined as 

5«(X1,X2) = (/»(X0/«(X2)>, (5) 

is one of the most important statistical descriptors of random media. It also can be inter- 
preted as the probability that two randomly chosen points xi and X2 both lie in phase i. 
For statistical homogeneous and isotropic media, only depends on scalar distances, i.e., 

^«(xi,X2) = 5«(|r|), (6) 

where r = X12. 

Global information about the surface of the ith phase may be obtained by ensemble 
averaging the gradient of J^*^(x). Since V/*-*^(x) is different from zero only on the interfaces 
of the ith phase, the corresponding specific surface defined as the total area of the interfaces 
divided by the volume of the medium is given by [l| 

.. = (|V/»(x)|>. (7) 

Note that there are other higher-order surface correlation functions which are discussed in 
detail by Torquato 

The calculation of higher-order correlation functions encounters both analytical and nu- 
merical difficulties, and very few experimental results needed for comparison purposes are 
available so far. However, their importance in the description of collective phenomena is in- 
disputable. A possible pragmatic approach is to study more complex lower-order correlation 
functions; for instance, the two-point cluster function C*^*-'(xi, X2) defined as the pro bability 
that two randomly chosen points xi and X2 belong to the same cluster of phase i [SOt]; or the 
lineal-path function L^*^(xi,X2) defined as the probability that the entire line segment be- 



tween points xi and X2 lies in phase i (311 • C-^H^i) ^2) and L'^*)(xi,X2) of the reconstructed 



media are sometimes computed to study the non-uniqueness issue of the reconstruction 
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III. NECESSARY CONDITIONS ON THE TWO-POINT CORRELATION FUNC- 
TION 



The task of determining the necessary and sufficient conditions that (r) must possess 
is very complex. In the context of stochastic processes in time (one-dimensional processes), 
it has been shown that the autocovariance functions must not only meet all the necessary 
conditions we will present in this section but another condition on "corner-positive" matrices 



32l |. Since little is known about corner-positive matrices, this theorem is very difficult to 



apply in practice. Thus, when determining whether a hypothetical function is realizable 
or not, we will ffist check all the necessary conditions collected here and then use the 
construction technique to generate realizations of the random medium associated with the 
hypothetical function as further verification. 



19 



20|. For a two- 



A. Known Necessary Conditions 

Here we collect all of the known necessary conditions on 5*2 
phase statistically homogeneous medium, the two-point correlation function for phase 2 is 
simply related to the corresponding function for phase 1 via the expression 

^f(r) = S«(r)- 201 + 1, (8) 
and the autocovariance function 

X(r) = S«(r)- 01^ = 5?) (r)-02^ (9) 
for phase 1 is equal to that for phase 2. Generally, for r = 0, 

5«(O) = 0„ (10) 
and in the absence of any long-range order. 



hm 5«(r)-.0,l (11) 

r +00 



An important necessary condition of realizable 5*2*'' (r) for a two-phase statistically ho- 
mogeneous medium with d dimensions is that the d-dimensional Fourier transform of the 
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autocovariance function x(r), denoted by x(k) must be non-negative for all wave vectors 
i.e., for all k 



X(k) = / x(r)e-^''-Mr > 0. (12) 

This non-negativity result is sometimes called the Wiener-Khint chine condition, which phys- 
ically results since x(k) is proportional to the scattered radiation intensity. The two-point 
correlation function must satisfy the following bounds for all r 

O<S«(r)<0,, (13) 
and the corresponding bounds on the autocovariance function are given by 

- min(0i^ 02') < X(r) < 0102- (14) 

A corollary of Eq. iHM recently derived by Torquato 13] states that the infimum of any 
two-point correlation function of a statistically homogeneous medium must satisfy the in- 
equalities 

max (0, 20i - 1) < inf \s^^ (r)j < 0^. (15) 

Another necessary condition on 5*2*^ (r) in the case of statistically homogeneous and 
isotropic media, i.e., when S'2*''(r) is dependent only the distance r = |r|, is that its derivative 
at r = is strictly negative for all < 0j < 1: 

^l^=»^^l..=o<0. (16, 

This is a consequence of the fact that slope at r = is proportional to the negative of 
the specific surface Taking that it is axiomatic that S'2*''(|r|) is an even function, i.e., 
'S'2*'*(kl) = *S'2*''(— |r|), then it is non- analytic at the origin. 



A lesser-known necessary condition for statistically 
"triangular inequality" that was first derived by Shepp 
eron jssll: 



lomogeneous media is the so-called 
and later rediscovered by Math- 



S^"{r)>S'i>{s) + S^'>{t)~<P,, (17) 
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where r = t — s. Note that if the autocovariance x(r) of a statistically homogeneous and 
isotropic medium is monotonically decreasing, nonnegative and convex (i.e., d^x/d^?" > 0), 
then it satisfies the triangular inequality Eq. (I17p . The triangular inequality implies several 
point-wise conditions on the two-point correlation function. For example, for statistically 
homogeneous and isotropic media, the triangular inequality implies the condition given by 
Eq. (fT6!) . the fact that the steepest descent of the two-point correlation function occurs at 
the origin 2o| : 



d5«(r) 



dr 



\r=0 



> 



dr 



(18) 



and the fact that 5*2* ""(r) must be convex at the origin |34l |: 



2 nii) 



d's: 



2 , _ d'x, 

r=0 — -r-^\r=0 ^ U. 



dr^ " dr^ " ^ 
Torquato [19] showed that the triangular inequality is actually a special case of the more 
general condition : 



m m 



i=l j=l 



(20) 



where Ei = ±1 (z = l,...,rri and m is odd). Note that by choosing m = 3; 6162 = 1, 
EiEs = £2^3 = —1, Eq- (fl^ can be rediscovered, li m = 3; 6162 = Sie^ = £2^3 = 1 are chosen 
instead, another "triangular inequality" can be obtained, i.e.. 



S^^ (r) > -S^^ (s) - S^^ (t) + - 0,) , (21) 

where r = t — s. Equation (pT!) was first derived by Quintanilla [ssl. 

Equation fl20l) is a much stronger necessary condition that implies that there are other 
necessary conditions beyond those identified thus far. However, Eq. (!20!) is difficult to check 
in practice, because it does not have a simple spectral analog. One possible method is to 
randomly generate a set of m points and compute the value of Xij = xi^i^'^^j)- Among these 
values of Xij, select the largest m ones and set their coefficients SiEj equal to —1. Thus, we 
have m equations for m e^'s. Then we can substitute the solved e^'s into Eq. fl20l) and check 
the inequality. If the inequality holds, then we can generate several different sets of random 
points and test the inequality in the same way. 



10 





FIG. 1: Numerical support of the conjecture: (a) Two-point correlation functions of black phase 

for hypothetical and constructed medium, (b) Constructed medium for which S2 best matches 
the target one. The linear size of the system N = 200 (pixels), volume fraction of black pixels 
4)1 = 0.227. 

B. Conjecture On A New Necessary Condition 



Besides the aforementioned explicit necessary conditions, we find in two-dimensional sim- 
ulations that the value of the second peak of a nonmonotonic S2{r) for a statistically homo- 
geneous and isotropic two-dimensional medium is always smaller than that of the medium 
composed of circular disks on a triangular lattice at fixed volume fractions, i.e.. 
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Si'\r,,4>,)<S^\r'^-^r\ (22) 

where Vp and denote the positions of the second peak for the two media, respectively 
and the superscript "tri" denotes the medium composed of disks on a triangular lattice. 
A hypothetical damped-oscillating two-point correlation function with an artificially higher 
second peak than that of the medium composed of disks on a triangular lattice at fixed 
volume fractions is tested by the construction algorithm. The results in Fig. [1] show that 
the structure for which S*^ ( "i?" denotes the black phase) best matches the target function 
indeed has its "particles" arranged on triangular lattice while the second peak of the target 
still cannot be reached. 

Here we make the conjecture that for any d-dimensional statistically homogeneous and 
isotropic medium with a two-point correlation function 5*2 (r) that is nonmonotonic in r, the 
value of the first peak of its 5*2*'' (r) away from origin is bounded from above by the value of 
the first peak of the two-point correlation function associated with the densest packings of 
(i-dimensional identical hard spheres at fixed a volume fraction, i.e.. 



5?^(r,;0.)<4^Hr;;0,)'^^^ (23) 

where the superscript "cps" denotes closest packings of spheres, for the first three dimensions, 
they are regular array of hard rods, hard disks on triangular lattice, and hard spheres on 
face-centered cubic lattice, respectively. For = 4 and d = 5, the densest packings are 
believed to be four- and five-dimensional checkerboard lattice packing, respectively [36|. 
Note that this conjectured condition could be a corollary of Eq. (1201) or some other unknown 
necessary conditions. 



IV. MODELING TWO-POINT CORRELATION FUNCTION VIA BASIS FUNC- 
TIONS 

A. Combination of Realizable Two-Point Correlation Functions 

It is first shown by Shepp [2^ that the convex combination and product of two realizable 
scaled autocovariance functions for one- dimensional statistically homogeneous media satisfy 
all known necessary conditions, i.e.. 
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f%r) = ai/i(r) + a2/2(r) 



(24) 



/nO = /iW/2(r), 

where < < 1 (i = 1, 2), ai + ^2 = 1 and the superscripts "c" and "p" denote "combina- 
tion" and "product", respectively. The scaled autocovariance function /(r) of a statistically 
homogeneous material is defined as 19| 



X[r) 



5«(r) 



(25) 



The necessary conditions for realizable scaled autocovariance function /(r) can be ob- 
tained from Eq. (!25|) and the equations by which the necessary conditions for realizable 
two-point correlation function 5*2*'' (r) are given. From Eqs. ( ITTI) and ( 12T|) . we can obtain the 
triangular inequalities for /(r), respectively, 



/(r)>/(s) + /(t)-l. 



(26) 



/(r) > -/(s) - /(t) - 1. 
Moreover, the bounds of /(r) become 



(27) 



— mm 



< /(r) < 1, 



(28) 



and the corollary Eq. f|T5l) is equivalent to 



— mm 



1 ll 

'^2 01 



< /inf < 0, 



(29) 



where /inf is the infimum of /(r). Our focus in this paper will be hypothetical continuous 
functions /(r) that are dependent only on the scalar distance r = |r| which could potentially 
correspond to statistically homogeneous and isotropic media without long range order, i.e.. 



/(O) = 1, hm /(r) 0. (30) 

r^oo 

/(r) is also absolutely integrable so that the Fourier transform of /(r) exists and is given by 
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oo 

f{k) = i2.r>^l/-'f(r):!m^dr > 0, (31) 



where A; = |k| and J„(x) is the Bessel function of order v. 

Generahzation of Eq. to higher dimensions is straightforward. Suppose fi{r) {i = 
1, m) are the scaled autocovariance functions for (i-dimensional statistically homogeneous 
and isotropic media, then the convex combination /^(r) and product /^(r) defined as 

m 
i=l 

(32) 

m 

F{r) = n Mr), 

i=l 

satisfy all known necessary conditions, where < ai < 1 {i = l,...,m) and Yl^i^i — ^■ 
Equation fl32l) is of great fundamental and practical importance. On the one hand, it enables 
us to construct new realizable two-point correlation functions with properties of interest, cor- 
responding to structures of interest, from a set of known functions. Thus, one can categorize 
microstructures with the set of known functions and the proper combinations. On the other 
hand, suppose that we can find a "full" set of those basis scaled autocovariance functions 
{fi{r)}iLi, then the scaled autocovariance function of any statistically homogeneous and 
isotropic medium can be expressed in term of the combinations of the basis functions, i.e., 

fir) = p[{/.(r)}r=i] = p[/i(r), /2(r), /^(r)], (33) 

where p denotes a map composed of convex combinations and products of {/i('^)}iILi- For 
example, for m = 5, a possible explicit form for p is 

P[{Mr)}U = (r) + «2[/5i/2(r) + /?2/3(r)] + as[h{r)Mr)], (34) 

where < ai, Pj < 1 and '^^di = Ylijl^j = 1 = l,2,3;j = 1,2). Once the scaled 
autocovariance function (or equivalently the two-point correlation function) of a medium is 
known, an effective reconstruction procedure enables one to generate accurate structures at 
will, and subsequent analysis can be performed on the image to obtain desired macroscopic 
properties of the medium. In other words, the medium is actually modeled by a set of 
basis scaled autocovariance functions {fiir)}'^^ and a particular map p[{fi{r)}'^^. There 
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could be different choices of the basis functions (hke different basis choices of a Hilbert 
space), and we would like the basis functions to have nice mathematical properties, such 
as simple analytical forms. Let {fi{r)}'^i denotes our choice of the basis functions. Thus, 
the media can be represented merely by different maps p^'s. Note that a hypothetical two- 
point correlation function corresponds to a hypothetical map p° effective construction 
algorithms can be used to test the realizability of 



B. Choice of Basis Functions 

A systematic way of determining the basis functions {ff}]li is not available yet. Here we 
take the first step to determine the bases by considering certain known realizable analytical 
two-point correlation functions and the corresponding scaled autocovariance functions. For 
convenience, we categorize these functions into three families: {i) monotonically decreasing 
functions; {ii) damped- oscillating functions; and {Hi) functions of known constructions. 

The family of monotonically decreasing functions includes the simple exponentially de- 



creasing function introduced by Debye 17| and polynomial functions. The former is given 
by 



= exp(-r/a), r > 0, 



(35) 



where a is a correlation length, corresponding to structures in which one phase consists of 



"random shapes and sizes" 
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37( 1 (shown in Fig. [2]). It is now known that certain types of 



space tessellations have autocovariance functions given by Eq. (!35l) [381]. We have referred 



to this class of structures as Debye random media 

Another example of monotonically decreasing functions is the family of polynomials of 
order n {n > 1) given by 



fpir) 



(1-r/a)" 0<r<a, 



(36) 



r > a, 

where a is the correlation length. The polynomial function of order 1 is shown to be realizable 
only for a statistically homogeneous two-phase medium in one dimension [l^, [20I. We 
have constructed for the first time realizations of random media in dimensions d < n that 
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FIG. 2: (a) Debye random medium function /^(r) with a = 5. (b) A realization of Debye random 
media with the volume fractions 0i = 0.68, 02 = 0.32. 

correspond to the polynomial function of order n with very high numerical precision using 
the Yeong-Torquato construction technique 2|]. However, for dimensions higher than n, 
Eq. ( !36l) violates Eq. ( l3TI) . We will henceforth assume that the polynomial functions of 
order n are realizable in dimensions d < n. 

An example of the family of damped-oscillating functions is given by 

foir) = ^^Ai&cp{-'r/ai)cos{qir + r > 0, (37) 



16 



where the parameters Ai and Oi {i = 1,2,...) control the amphtude of the fo profile, qi is the 
wavenumber and ipi is the phase angle. In general, the wavenumber should be a function of 
r which could correspond to the different distances between successive neighbor shells of a 
crystalline material. Note that qi and ipi need to be carefully chosen such that foi'f') satisfies 
all known necessary conditions. 





(b) 

FIG. 3: (a) Scaled autocovariance function fsir) of two-dimensional identical overlapping disks 
with volume fractions (pi = 0.45, 02 = 0.55. (b) A realization of two-dimensional identical over- 
lapping disks with volume fractions cpi = 0.45, (p2 = 0.55. The radius of disks R = 5. 



The family of functions of known constructions includes scaled autocovariance functions 
of (i-dimensional identical overlapping spheres |l|, Q and symmetric-cell materials |l|, Q . 
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For overlapping spheres of radius R, the scaled autocovariance function for the particle phase 
(spheres) is given by 



= eM-pv^(r.R)]-<Pl (38) 

0102 

where 0i and 02 are volume fractions of the spheres and matrix respectively, p = N/V is the 
number density of spheres, and V2{r; R) is the union volume of two spheres of radius R whose 
centers are separated by r. For the first three space dimensions, the latter is respectively 
given by 



vi{R) 



29 (r - 2R) + (l + ^) 0(2i? - r), (39) 



= 2e(r - 2i?) + 2 

Vi{R) ^ ' TT 



TT 



r 

2R 



1 - 



— cos 



\2R} 



e{2R-r), (40) 



V2ir;R) 
viiR) 



2e(r - 2R) + 



1 + 



3r 
4R 



r \ 3 

16 v^y 



e(2i?-r) 



(41) 



where G(x) is the Heaviside function, and vi{R) is the volume of a c?- dimensional sphere of 
radius R given by 



vi{R) 



TT 



d/2 



(42) 



r(i + rf/2)' 

where T{x) is the gamma function. For c? = 1,2 and 3, vi{R) = 2R,ttR^ and 4ttR^/3, 
respectively. 

Two-phase symmetric-cell materials are constructed by partitioning space into cells of 
arbitrary shapes and sizes, with cells being randomly designated as phase 1 and 2 with 
probability 0i and 02 from a uniform distribution [l| . For such a statistically homogeneous 
and isotropic medium, the scaled autocovariance function is given by 



fc{r)=wi'\r), (43) 

where W2^\r) is the probability that two points separated by distance r are in the same 
cell. This quantity is only a function of cell shapes and sizes, depending on a dimensionless 
size-averaged intersection volume of two cells defined by 
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(b) 



FIG. 4: (a) Scaled autocovariance function fcif) of two-dimensional random checkerboard, (b) A 
realization of two-dimensional random checkerboard. The length of the cells a = 10. 



{vf\r;R)) 



R 



{vi{r))R 

where R is the size parameter for each cell, {v^*{r; R))^^ is the size-averaged intersection 
volume of two cells whose centers are separated by r and is the size-averaged single- 

cell volume. The random checkerboard is a very useful model of symmetric-cell material 
because its fc{i^) is known analytically for d = 1 and 2 For ci = 1, it is easy to verify 
that the probability of finding two points in the same one-dimensional cell is given by 
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1 — r/a, < r < a 



(45) 



0, r > a. 

where a is the length of the side of a square cell. Note that Eq. is just the polynomial 
function of order one [cf. Eq. fl36l) ] for one-dimensional homogeneous media. For d = 2, 
VTg"^^ {^) is given by 



1 - 



TT 



1 + 


1 


/r 








TT 


\a 






4 




+ 








TT 



- -4 



< r < a, 



r \ 



0, 



— 1 — cos 
a/ 



> V2a. 



<r< V2a, (46) 



For d = 3, W2^\r) is given by 



W2^i)(r) = - / / Vr(r,^,0)sin 



vr 



^0 



where 



(47) 



Vr(r, ^, 0) = [1 — r cos [1 — r sin sin 0] [1 — r sin cos 0] 

(48) 

xO(l — rcos6')B(l — r sin^sin 0)0(1 — rsin6'cos0). 

Note that ^"^ i^) d = 3 does not have a simple analytical form. 

There are some other known scaled autocovariance functions for different types of materi- 

n 

als (e.g., (i-dimensional identical hard spheres [1|) and realizable hypothetical functions (e.g., 
complementary error function, see Appendix). However we do not include these functions 
in our basis function set because they do not have simple analytical mathematical forms. 
Note that some of the basis functions are dimension-dependent (e.g., fo and fc), and thus 
the proper forms of basis functions should be used for different dimensions. In the following 
discussion, without further specification, we will focus on two-dimensional cases. 
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V. GENERATING REALIZATIONS OF HETEROGENEOUS MATERIALS 



Consider a digitized (i.e., pixelized) representation of a heterogeneous material. Different 
colored pixels (in a discrete coloring scheme) may have numerous interpretations. The image 
can reflect different properties, such as the geometry captured by a photographic image, 
topology of temperature and scalar velocity fields in fluids, distribution of magnitudes of 
electric and magnetic fields in the medium, or variations in chemo-physical properties of the 
medium. In the last case, typical examples are composite materials in which the different 
phases may have different thermal, elastic or electromagnetic properties, to name a few. 

A. Exact Equations for Digitized Media 

Our focus in this paper is two-dimensional, two-phase statistically homogeneous and 
isotropic random media composed of black and white pixels. Such a system can be repre- 
sented two-dimensional array, i.e., 

hi h2 
J hi h2 

Ini Im 

where the integer categorizes the linear size of the system (A^^ is the total number of 
pixels in the system) and the entries hj {hi = 1, can only take the value of or 1, 

which correspond to the white and black phases, respectively. Note that Eq. ( H9|) is only an 
abstract representation and the real morphological configuration of the medium also depends 
on the choice of lattices. For example, as shown in Fig. [5l the isotropic medium composed 
of overlapping disks generated on a square lattice and the anisotropic medium composed 
of orientated ellipses generated on a triangular lattice have the same array presentation. A 
vector distance in the digitized medium can be uniquely expressed as 

r = niei + n2e2, (50) 

where ei and e2 are lattice vectors for the particular lattice and ni, n2 are integers. For 
example, for square lattice, ei = i, 62 = j, where i, j are unit vectors along horizontal and 



hN 
hN 

Inn 



(49) 
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vertical directions, respectively; while for triangular lattice ei = ^/3i + j/2, 62 = j. 




(a) 



(b) 



FIG. 5: Different digitized media witli tlie same arrary representation: (a) Overlapping disks 
generated on square lattice with a square unit cell, (b) Orientated overlapping ellipses generated 
on triangular lattice with a rhombical unit cell. 

Without loss of generahty, we choose the black phase to be the phase of interest and 
assume periodic boundary condition is applied, which is commonly used in computer simu- 
lations. The two-point correlation function 6*2 (r) of the black phase can be calculated based 
on its probabilistic nature, i.e., the probability of finding two points separated by the vector 
distance r in the same phase. The value of two-point correlation function for a particular 
r = niGi + 71,262 is given by 



where lij are entries of I defined in Eq. ( H9l) . ni and n2 are integers satisfying |ni|, |n2| < 
[A^/2] due to minimum image distance convention. For statistically isotropic media, the 
two-point correlation function only depends on the magnitude of r, i.e., r = |r|, thus we 
have 




(51) 




(52) 



where 
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n= {{m,n) \ + n'^ = r^, r <[N/2]} , (53) 

and u is the number of elements of set Q. 

It is well known that the two-point correlation function cannot completely specify a 
two-phase heterogeneous material alone. Here, we provide a proof of this statement for two- 
dimensional statistically homogeneous digitized media. Note that the proof trivially extends 
to any dimension. Suppose we already know the value of 5*2 (r) for every vector distance r, 
using Eq. flSTl) . we can obtain a set of equations of lij, i.e., for each r = riiBi + 72262 

N N 

5Z 5Z ^iihi+ni)(j+n2) " A^^5'2(ni, T^s) = 0, (54) 
i=l j=l 

Since the digitized medium is represented by I, once we obtain all the entries lij (the 
unknowns in Eq. fl5^ ). the medium is (re) constructed. In Eq. fl5^ . the number of unknowns 
Nu equals A^^ and the number of equations A^e equals the number of all possible vector 
distances in the digitized medium. To calculate Ng, all possible different combinations of 
integers rii and n2 subjected to |ni|, |?t,2| < [N/2] need to be considered. This quantity is 
given by Nf, = 4[A^/2]^ — 8[A^/2] -|-6, which is smaller than the number of unknowns = A^^ 
for normal-sized systems (i.e., = 10 ~ 10^). 

Similar proof is applied to the case when we average over the angles of vector r to yield 
5*2 (r) that depends only on the radius distance r. The angle-averaged equations of lij are 
given by 



' N N 
.i=l j=l 



ujN^S2ir) = 0, (55) 



(m,n)Gf! 

where Q is given by Eq. ( |53l) . The number of equations is even smaller for the angle-averaged 
case, while the number of unknowns does not change. The analysis shows that one could 
never find a unique solution of lij from neither Eq. (15^ nor Eq. ( l55ll unless some assumptions 
have been made that reduce Nu such that A^^ = N^. For example, in an interesting model to 
study spatial distribution of algae, the algae are put on top of each other in order to reduce 
the unknowns 4l(]. In general cases, one could use the stochastic optimization procedure 
(i.e., simulated annealing) to find solutions of Eq. ( !54l) and Eq. ( !55|) . Note that 

although the aforementioned proofs focus on two-dimensional media, they trivially extend 
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to any dimensions ( e.g., d = 1,3). 



B. Stochastic Optimization Procedure 

Generally, consider a given set of correlation functions /°(ri, r2, r„) of the phase of 
interest that provides partial information on the random medium. The index a is used to 
denote the type of correlation functions. Note that the set of should not be confused 
with the basis function set p, the former contains correlation functions of different type, i.e., 
two-point correlation function, lineal-path function, two-point cluster function, etc., while 
the latter contains basis functions through which the scaled autocovariance function of the 
medium of interest can be expressed. The information contained in /" could be obtained 
either from experiments or it could represent a hypothetical medium based on simple models. 
In both cases we would like to generate the underlying micro-structure with a specified set 
of correlation functions. In the former case, the formulated inverse problem is frequently 
referred to as a "reconstruction" procedure, and in the latter CclSG clS db "construction" . 



As we have noted earilier, it is natural 



to formulate the construction or reconstruction 



problem as an optimization problem . The discrepancies between the statistical 

properties of the best generated structure and the imposed ones is minimized. This can 
be readily achieved by introducing the "energy" function E defined as a sum of squared 
differences between target correlation functions, which we denote by and those calculated 
from generated structures, i.e., 

^= E E[/n(^i'^2,...,r.)-/:(ri,r2,...,r„)]'. (56) 

ri,r2,...,r„ a 

Note that for every generated structure (configuration), these is a set of corresponding If 
we consider every structure (configuration) as a "state" of the system, E can be considered 
as a function of the states. The optimization technique suitable for the problem at hand 
is the method of simulated annealing [igI. It is a popular method for the optimization 
of large-scale problems, especially those where a global minimum is hidden among many 
local extrema. The concept of finding the lowest energy state by simulated annealing is 
based on a well-known physical fact: If a system is heated to a high temperature T and 
then slowly cooled down to absolute zero, the system equilibrates to its ground state. At 
a given temperature T, the probability of being in a state with energy E is given by the 
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Boltzmann distribution P{E) ~ exp{—E/T). At each annealing step k, the system is allowed 
to evolve long enough to thermalize at T{k). The temperature is then lowered according to 
a prescribed annealing schedule T{k) until the energy of the system approaches it ground 
state value within an acceptable tolerance. It is important to keep the annealing rate slow 
enough in order to avoid trapping in some metastable states. 

In our problem, the discrete configuration space includes the states of all possible pixel 
allocations. Starting from a given state (current configuration), a new state (new configu- 
ration) can be obtained by interchanging two arbitrarily selected pixels of different phases. 
This simple evolving procedure preserves the volume fraction of all involved phases and 
guarantees ergodicity in the sense that each state is accessible from any other state by a 
finite number of interchange steps. However, in the later stage of the procedure, biased and 
more sophisticated interchange rules, i.e., surface optimization, could be used to improve the 
efficiency. We choose the Metropolis algorithm as the acceptance criterion: the acceptance 
probability P for the pixel interchange is given by 



where AE = Enew~Eoid- The temperature T is initially chosen so that the initial acceptance 
probability for a pixel interchange with AE > averages approximately 0.5. An inverse log- 
arithmic annealing schedule which decreases the temperature according to T{k) ~ l/ln(A;) 
would in principle evolve the system to its ground state. However, such a slow annealing 
schedule is difficult to achieve in practice. Thus, we will adopt the more popular and faster 
annealing schedule T{k)/T{0) = X'^, where constant A must be less than but close to one. 
This may yield suboptimal results, but, for practical purposes, will be sufficient. The con- 
vergence to an optimum is no longer guaranteed, and the system is likely to freeze in one of 
the local minima if the thermalization and annealing rate are not adequately chosen. 

The two-point correlation function of a statistically homogeneous and isotropic medium 
is the focus of this paper. In this case, Eq. fl56l) reduces to 




V 




1 



AE < 0, 



(57) 




(58) 
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S'2-sampling algorithm. One of the most commonly used and efficient S'2-sanip^ 
is the orthogonal-sampling algorithm introduced by Yeong and Torquato 



Since for every configuration (structure), the corresponding two-point correlation function 
needs to be computed, the efficiency of the construction or reconstruction is mainly de- 
termined by the efficiency of the S'2-sampling algorithm. Furthermore, the properties of 
generated configurations (structures), i.e., isotropy of the medium, is also affected by the 

ing algorithms 
Due to the 

isotropic nature of the medium, every sampling direction should be equivalent. For sim- 
plicity, two orthogonal directions (usually the horizontal and vertical directions of a square 
lattice) are chosen and the two-point correlation function is sampled along these directions 
and averaged. At each pixel interchange, only the values of 5*2 (r) sampled along the rows 
and columns that contain the interchange pixels are changed. Thus, the complexity of the 
algorithm is 0{N), where is the linear size of the system. However, for certain media 
with long-range correlations, the generated media have microstructures with two orthogonal 
anisotropic directions due to the biased sampling. Modifications of the orthogonal-sampling 
algorithm to preserve the isotropy of the underlying medium have been proposed, such as 
adding more sampling directions and using more isotropic lattices {4, 5|. Cule and Torquato 
i ntroduced a new isotropy-preserving fast Fourier transform (FFT) algorithm [4]. At each 
pixel interchange step, the two-point correlation function 5*2 (r) containing angle informa- 
tion is calculated in momentum space using an efficient FFT algorithm. Since information 
of all directions is considered, the generated media always have the required isotropy struc- 
tures. However, since the complexity of FFT is 0(A^log2 A^), the algorithm is relatively time 
consuming. 

We have developed an efficient and isotropy-preserving algorithm, namely, the Lattice- 
Point algorithm by considering the black pixels as hard "particles" on a particular lattice. 
The two-point correlation function is then computed in a similar way of obtaining the pair 
correlation function g2{r) for an isotropic point process At each Monte-Carlo step, the 
randomly selected "particle" (black pixel) is given a random displacement subjected to the 
nonoverlapping constraint and the distances between the moved "particle" and all the other 
"particles" need to be recomputed. Thus the complexity of the algorithm is 0{N). Since 
all directions are effectively sampled, constructions based on the angle-averaged 5*2 (r) well 
preserve isotropy of the media. A detailed discussion and applications of the algorithm will 
be included in the second paper of this series 29|]. In this paper, we only provide several 
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illustrative examples generated from this algorithm. 

VI. ILLUSTRATIVE EXAMPLES 

As illustrative examples, we use the aforementioned (re) construction techniques to in- 
vestigate both deterministic, crystal-like structures, and random systems. We also study 
a hypothetical medium with the two-point correlation function obtained from convex com- 
bination of known ones. In the case of completely deterministic structures, the algorithm 
produces almost perfect reconstructions. However, the optimization of disordered struc- 
tures is significantly harder. Furthermore, we will see that for the media with long-range 
correlations, e.g., a damped-oscillating 5*2 (r), the orthogonal-sampling algorithm may pro- 
duce unexpected anisotropy, while the Lattice-Point algorithm well preserves isotropy of the 
media. 

A. Regular Array of Nonoverlapping Disks 




FIG. 6: A realization of square array of nonoverlapping disks. The linear size of the system = 200 
pixels, volume fraction of black pixels (pi = 0.326. 



First, we consider specific two-dimensional and two-phase structure composed of a square 
array of nonoverlapping disks, as shown in Fig. El This morphology may be viewed as a 
cross section of two-phase materials containing rod- or fiber-like inclusions. Various trans- 
port properties of these materials have been well explored because of their practical and 



theoretical importance in materials science 



42|. 
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The regular structure is discretized by introducing an x square lattice. The volume 
fractions of black and white phases are 0i = 0.326, 02 = 0.674, respectively. The target 
two-point correlation of the digitized medium is sampled using both the orthogonal and the 
Lattice-Point algorithm for comparison purpose. The simulations start from random initial 
configurations (i.e., random collections of black and white pixels), at some initial temperature 
To, with fixed volume fractions 0j. At each Monte-Carlo (MC) step, when an attempt 
to exchange two randomly chosen pixels with different colors (or to randomly displace a 
chosen black pixel) is made, 5*2 (r) is efficiently recomputed by using the orthogonal-sampling 
algorithm (or the Lattice- Point algorithm). The set of constants {Xmc, ^tot, A} specifies the 
annealing schedule: At each temperature, the system is thermalized until either Xmc^"^ MC 
moves are accepted or the total number of attempts to change the original configurations 
reaches the value XtotN"^- Subsequently, the system temperature is decreased by the reduction 
factor A, i.e., T^ew = XTom- 




(a) (b) 

FIG. 7: Reconstructed structures: (a) Square array of almost circular particles generated by 
the Lattice-Point algorithm, (b) Square array of particles generated by the orthogonal-sampling 
algorithm. 

The reconstruction results are shown in Fig. [71 Both of the algorithms are able to re- 
produce the exact global square-array arrangement of clusters of black pixels. This implies 
that the two-point correlation function of regular configurations contains enough structural 
information to properly characterize the long-range correlations. However, it is clear that 
the structure generated by the Lattice-Point algorithm has a better local arrangement of the 
pixels (i.e., the shape of the particles) than that generated by the orthogonal-sampling algo- 
rithm. This is because the orthogonal algorithm only uses structural information along two 
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directions, which is not sufficient to reproduce detailed local structures, while the Lattice- 
Point algorithm efficiently uses information along all possible directions. 

B. Hypothetical Random Media with Long-Range Correlations 

In this example, we will generate two-dimensional statistically homogeneous and isotropic 
random media with long-range correlations (i.e., nontrivial inter-particle interactions). Ex- 
amples of this type of media include low-density fluids and amorphous materials (i.e., porous 
media, randomly polymerized plastics, glass, etcj. A meaningful, yet nontrivial, two-point 
correlation function capturing these features is 0, [s], Q 

52(r) = 0? + 0i02e-^/'-o^^, (59) 

kr 

where k = 27r/ao. Here tq and oq are two characteristic length scales. The overall exponential 
damping is controlled by the correlation length ro, determining the maximum correlations 
in the system. The constant ag determines oscillations in the term sm{kr)/{kr) which also 
decays with increasing r, such that ao can reduce the effective range of tq. Interestingly, this 
hypothetical function is not exactly realizable, because it violates the convexity condition 
Eq. f|T9|) at the origin, or more generally the triangular inequality Eq. f|T7|) . However, we 
mainly focus on its damped-oscillating property and we will see that the algorithms are 
robust enough to detect violation of the convexity condition. 

For comparison purposes, both the orthogonal-sampling algorithm and the Lattice-Point 
algorithm are used in the construction, the results are shown in Fig. [HI At a lower density of 
the black phase 0i, ao is manifested as a characteristic repulsion among different elements 
with diameter of order oq. The repulsion vanishes beyond the length scale tq. At a higher 
density, both length scales oq and ro are clearly noticeable in the distribution of the black 
and white phases. Note that the structures generated by the orthogonal-sampling algorithm 
exhibit some anisotropy features, i.e., containing stripes along ±45 degree directions, which 
implies that the orthogonal-sampling algorithm should be used with care in the case where 
the medium has long-range correlations. 

The target two-point correlation function 5*2 (r) for (pi = 0.2 and 5*2 (r) sampled from 
generated structures are shown in Fig. [9l It can be seen clearly that 5*2 (r) is non-convex 
at the origin and the largest discrepancies between 5*2 (r) and 5*2 occur around the origin 
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FIG. 8: (a) Media with 5*2 (r) given by Eq. (j59p generated by the orthogonal-samphng algorithm: 
Left panel, volume fraction of black pixels (pi = 0.2. Right panel, volume fraction of black pixels 
01 = 0.5. The linear size of the systems N = 200. (b) Media with S2{r) given by Eq. (|59p generated 
by the Lattice-Point algorithm: Left panel, volume fraction of black pixels (pi = 0.2. Right panel, 
volume fraction of black pixels (pi = 0.5. The linear size of the systems N = 200. 

because 5*2 (r) satisfies Eq. f|T9|) . This implies that our algorithms are robust enough and 
can be used to test realizability of hypothetical functions. Note that the following classes of 
functions UM 



f{r) = exp 



^ \ Of 

a 



a > 1, 



(60) 



and 
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FIG. 9: Target two-point correlation function given by Eq. (j59p and that of constructed media 
with volume fraction (pi = 0.2. 



cannot correpond to a two-phase medium in d dimensions also because of violation of trian- 
gular inequality Eq. 

C. Hypothetical Random Media with Realizable Correlation Functions 

In this last example, we study a hypothetical statistically homogeneous and isotropic 
medium whose scaled autocovariance function is a convex combination of Debye random 
medium function foir) and damped-oscillating function fo{^), i-e., 



where ai + a2 = 1. Since both fD{r) and fo{f) are independent of volume fractions, the 
medium with scaled autocovariance function /(r) has phase-inversion symmetry i.e., the 
structures with volume fraction of black pixels 0i = 0.2 is statistically the same with those 
having volume fraction of white pixels 02 = 0.2, if the colors of the two phases are inverted 
in the latter. 

Different constant pairs (ai, 0:2) can be used to construct /(r) with required properties. 




(61) 




(62) 
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FIG. 10: Combined scaled autocovariance function /(r) with coefficients ai = 0.25, 02 = 0.75. 




FIG. 11: Constructed media with scaled autocovariance function shown in Fig. [TUl (a) Volume 
fraction of black pixels 0i = 0.1. (b) Volume fraction of black pixels (pi = 0.3. (c) Volume fraction 
of black pixels 0i = 0.5. The linear size of the systems = 200. 

In particular, we choose two pairs: (0.25, 0.75) and (0.75, 0.25). The construction results 
obtained by application of the Lattice-Point algorithm are shown in Figs. [TOl [11] and Figs. \T2\ 
[T3I For tti = 0.25, a2 = 0.75, fo{f) is dominant in the combination. At lower densities, the 
generated structures resemble those with "pure" damped oscillationg two-point functions, 
i.e., dispersions of particles, altough they contain more clusters. At higher densities, some 
stripe-like structures and several (almost equal-sized) clusters can be identified. For ai = 
0.75, 02 = 0.25, foi^) is dominant in the combined f{r). Clusters with considerable sizes 
form even at low densities, which is a consequence of the large effective correlation length 
in foir). However, no stripe-like structures can be identified in the generated structures. 



32 



since the contribution of foi^) is significantly supressed by its damping nature and tlie small 
combination constants. 




50 



r (pixels) 



FIG. 12: Combined scaled autocovariance function /(r) with coefficients ai = 0.75, a2 = 0.25. 




FIG. 13: Constructed media with scaled autocovariance function shown in Fig. [H[ (a) Volume 
fraction of black pixels 0i = 0.1. (b) Volume fraction of black pixels (pi = 0.3. (c) Volume fraction 
of black pixels (pi = 0.5. The linear size of the systems N = 200. 

The results imply that even a simple combination of two basis functions enables one to 
obtain scaled autocovariance functions with properties of interest and to generate a variety 
of structures with controllable morphological features, e.g., local "particle" shape and cluster 
size. 
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VII. CONCLUSIONS 



In this paper, we have provided a general rigorous scheme to model and categorize two- 
phase statistically homogeneous and isotropic media. In particular, given a set of basis 
functions, we have shown that the medium can be modeled by a map p composed of convex 
combination and product operations. The basis functions should be realizable but, if they 
are not, they should at least satisfy all the known necessary conditions for a realizable 
autocovariance function. We have gathered all the known necessary conditions and made 
a conjecture on a possible new condition based on simulation results. A systematic way of 
determining basis functions is not available yet. We proposed a set of basis functions with 
simple analytical forms that capture salient microstructural features of two-phase random 
media. 

We give for the first time a rigorous mathematical formulation of the (re) construction 
problem and showed that the two-point correlation function alone cannot completely spec- 
ify a two-phase heterogeneous material. Moreover, we devised an efficient and isotropy- 
preserving (re) construction algorithm, namely, the Lattice-Point algorithm to generate real- 
izations of materials based on the Yeong-Torquato technique. We also provided an example 
of non-realizable yet non-trivial two-point correlation function and showed that our algo- 
rithm can be used to test realizability of hypothetical functions. An example of generating 
hypothetical random media with combined realizable correlation functions was given as an 
application of our general scheme. We showed that even a simple combination of two basis 
functions enables one to produce media with a variety of microstructures of interest and 
therefore a means of categorizing microstructures. 

We are investigating applications of our general scheme in order to model real materials. 
We are also developing more efficient (re) construction algorithms. There is a need for a 
theoretical and numerical analysis of the energy threshold of the algorithm, which is the 
aforementioned "acceptable tolerance". This quantity provides an indication of the extent 
to which the algorithms have reproduced the target structure and it is directly related to 
the non- uniqueness issue of reconstructions ^, ^,0,151. Moreover, additional realizable basis 
functions are needed to construct a complete basis set. Such work will be reported in our 
future publications. 
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APPENDIX A 



The complementary error function fce{f) is defined as 

9 /-oo 

Ur) = — / e^'^dt, (Al) 

Jr/a 

where r > and a is the effective correlation length. It is easy to check that fce{f) satisfies 
the known necessary conditions collected in this paper except for Eq. ( l20l) . which can only 
be checked for a finite number of cases. Realizations of the random medium associated 
with fcei'f') have been constructed using the Yeong-Torquato technique with very high nu- 
merical precision. Thus, we believe fce{f) to be a valid candidate for a realizable scaled 
autocovariance function. 
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